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Abstract. The band bending (BB) effect on the surface of the second-generation 
topological insulators implies a serious challenge to design transport devices. The BB 
is triggered by the effective electric field generated by charged impurities close to the 
surface and by the inhomogeneous charge distribution of the occupied surface states. 
Our self-consistent calculations in the Korringa-Kohn-Rostoker framework showed that 
in contrast to the bulk bands, the spectrum of the surface states is not bent at the 
surface. In turn, it is possible to tune the energy level of the Dirac point via the 
deposited surface dopants. In addition, the electrostatic modifications induced by the 
charged impurities on the surface induce long range oscillations in the charge density. 
For dopants located beneath the surface, however, these oscillations become highly 
suppressed. Our findings are in good agreement with recent experiments, however, our 
results indicate that the concentration of the surface doping cannot be estimated from 
the energy shift of the Dirac cone within the scope of the effective continuous model 
for the protected surface states. 


PACS numbers: 73.20.At; 

1. Introduction 

The theoretical discovery of the second-generation topological insulators [T] (2GTIs) 
triggered an intensive experimental effort to observe the predicted surface states 
[3l[2l|5llH[6l|7llHl|9l Ho] (SSs) being protected by time-reversal symmetry |TT|. It 
turned out, that the physical properties of the prepared samples are greatly affected 
by the electron acceptor/donor impurities, that can be found either in the bulk or on 
the surface. The created Bi 2 Se 3 samples are typically electron doped by the inner point 
defects [niHSlHlI. However, it has been shown that the Fermi level can be tuned by 
the insertion of further bulk dopants into the system [giiHi]. The presence of the 
charged impurities in the system and the inhomogeneous charge distribution close to 
the surface generates an effective electrostatic field, which can be probed experimentally, 
for example by second-harmonic generation [miiEi. 
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In addition, the evolved electrostatic field induces a band bending (BB) in the bulk 
band structure close to the surface, which was successfully observed by angle-resolved 
photomeisson (ARPES) experiments as well P El @1 EHl EZl E]. The experimental 
manipulation of the BB held was also accomplished by the insertion of bulk Cu dopants 
into the Bi 2 Se 3 matrix |9]. Moreover, ARPES experiments also demonstrated the 
possibility to shift the Dirac cone by applying charge dopants on the Bi 2 Se 3 surface 
[a E]. An experimental evidence for a large shift of the Dirac cone towards the 
conduction band was also reported by gated terahertz cyclotron resonance measurements 
performed on thin Bi 2 Se 3 hlm [T9]. 

Besides these comprehensive experimental studies, numerous theoretical works were 
also devoted to the description of the physical properties of the 2GTIs, including hrst 
principle calculations jll EQI EH EH EH EH EH]) tight binding [23l EH] El] or effective 
continuous m models. The structure of the 2GTIs can be described by a sequence of 
weakly bound quintuple layers (QLs), each consisting of hve atomic layers. The effect of 
the bulk dopants on the electronic behavior of Bi 2 Se 3 crystal was also addressed in recent 
theoretical studies |2H EH EH El] . In addition, Galanakis et. al. [2^ also suggested that 
the BB prohle and the energy of the Dirac point can be controlled by electrostatic effects. 
The importance of the BB has been demonstrated in other 2GTI based nanostructures 
as well, including topological/normal insulator interfaces [29l [SO] [311 122 ] • 

Motivated by these research Endings in this work we study the BB efiect and 
the Dirac cone shift on the surface of 2GTIs theoretically. In particular, we perform 
screened Korringa-Kohn-Rostoker (SKKR) calculations to examine the role of the 
charged dopants at Bi 2 Se 3 surface. The slow dynamics of the band bending process 
suggests that the charge accumulation at the surface is coupled to a much slower surface 
lattice relaxation |2] . In this work we do not aim to describe the time dependence of the 
outlined process, but rather to examine the efiect of the surface dopants on the band 
structure. In addition, according to Ref. 121 the lattice relaxation in the presence of 
adatoms is expected to be negligible small. Thus, in our surface calculations we consider 
a rigid lattice excluding any structural relaxation processes of the surface layers. 

The rest of the paper is organized as follows. In Sec. El we present the details of our 
numerical approach to study the bulk and surface properties of Bi 2 Se 3 . Than in Sec. [3] 
we examine the BB and the efiect of the surface dopants on the dispersion of the SSs. 
Finally we summarize our work in Sec. IH 

2. Details of the numerical calculations 

In order to describe the electron structure of Bi 2 Se 3 , we used the relativistic spin- 

polarized SKKR method [33]. The first-principles calculations were performed by 

density functional theory using the local spin-density approximation and the Geperley- 

Alder parametrization of the exchange correlation functional [3l] within the the atomic- 

sphere approximation (ASA). In our calculations we used an angular momentum cutoff 

/ = 2 
^max 
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The 2GTIs posses a rhombohedral lattice structure were the atoms are located in 
parallel layers forming a triangular lattice [3S].This lattice structure can be described 
by a periodic sequence of QLs, whereas each of the QL consists of hve strongly bound 
atomic layers. The QLs, on the other hand, are weakly bound to each other by van der 
Waals couplings. In particular, the QLs in the Bi 2 Se 3 crystal consist of atomic layers 
Sel-Bi-Se2-Bi-Sel, whereas Sel and Se2 are selenium atoms at inequivalent geometrical 
positions. Fig. [^a) shows the structure of one QL in the lattice. Since the QLs are 



Figure 1. (a) The structure of the Bi 2 Sea crystal within one QL. (b) The two- 

dimensional Brillouin zone related to the surface terminated by an atomic plane of 
a QL.(c) The scheme of the studied system. Six QLs are surrounded by a semi-infinite 
bulk of Bi 2 Se 3 crystal from the left and by vacuum from the right. 


weakly bound to each other, the crystal surface is favored to be terminated by Se atoms. 
Thus, in our calculations we assumed a flat surface formed by the last Se atomic plane 
of a QL. The considered surface has a hexagonal structure, with a 2D lattice constant 
a = 4.138 A and lattice vectors given by ai = (a, 0, 0) and a 2 = (— 0), where the 
axis z is perpendicular to the surface. The shortest period in the lattice structure along 
the z axis is given by three successive QLs PQ. Still, the periodicity of the Bi 2 Se 3 lattice 
in the 2 ; direction can be well described by a skew lattice vector a 3 = (5a, ^a, zql) 
(with Zql = 9.547 A standing for the hight of one QL) resulting in a period length of 
one QL only in the direction spanned by as. 

According to the SKKR method [33], we modeled the bulk system by a single 
QL surrounded by semi-inhnite bulk regions from both the left and right sides. The 
surface Greens functions of the semi-inhnite regions were calculated by an iterative 
method [36l[37|. The charge density distribution was calculated by means of the energy- 
dependent Greens function, integrated up to the Fermi energy. The position of the Fermi 
energy was determined to satisfy the total charge neutrality condition. In particular, 
because of the ASA used in our calculations, however, the value of the Fermi level is 
obtained by a numerical error. Though, the correct position of the Fermi energy is 
essential for insulator systems. Several previous works [38l [39l HQ] proposed a procedure 
to correct the value of the Fermi energy in a self-consistent way based on the Lloyd’s 
formula [TT]. Following this procedure we obtained the proper Fermi energy by re¬ 
normalizing the wave functions in order to obtain the correct space-integrated charge 
distribution. 






















Band bending at the surface of Bi 2 Se^ studied from first principles 


4 


For the surface calculations we considered an interface region surrounded by a semi¬ 
infinite bulk system from the left and by a vacuum from the right [see Fig. [11(c)]. The 
interface region was constructed from six QLs to support a smooth transition of the 
atomic potentials as we proceed from the bulk layers to the vacuum. We also included 
additional (in total eight) vacuum layers between the surface of the interface region and 
the semi-inhnite vacuum side. From an electrostatic point of view, the atomic potentials 
(as well as the charge density distribution) were obtained by imposing zero electrostatic 
held far away from the surface of the crystal. 

The described numerical method is sufficient to obtain plausible results for the 
studied system [22], however, the inclusion of empty spheres between the atomic 
layers further stabilized our numerical approach. In our calculations we used identical 


Type 

x|Al 

y[Al 

z[A] 

Rws [^] 

El 

0 

0 

0 

1.447 

Sel 

2.069 

1.195 

1.184 

1.587 

E2 

4.138 

2.389 

2.1408 

1.226 

Bi 

6.207 

3.584 

2.893 

1.758 

E3 

8.276 

4.778 

3.664 

1.1854 

Se2 

10.3450 

5.973 

4.773 

1.698 

E3 

12.414 

7.167 

5.883 

1.185 

Bi 

14.483 

8.362 

6.654 

1.758 

E2 

16.552 

9.556 

7.406 

1.226 

Sel 

18.621 

10.751 

8.363 

1.587 


Table 1. The x, y and z coordinates and the Wigner-Seitz radii {Rws) of the atomic 
and empty spheres in one QL. Notations Set and Se2 [El, E2, E3] stand for the 
inequivalent selenium [empty] spheres. The position of the other spheres in the lattice 
can be computed via the lattice vectors ai, sl 2 and as (see the text for details). 


geometrical parameters for the lattice structure as in Ref. [35], however, we optimized 
the positions and the Wigner-Seitz radii of the empty spheres to reproduce the main 
features of the experimentally observed band structure of the SSs.jl] Focusing on the 
direct band gap at the F point and on the slope of the Dirac cone, the obtained numerical 
parameters are summarized in Table [H 

Finally, the band structure can be obtained from the Bloch spectral function (BSF). 
For surface calculations the layer-resolved BSF depends on the energy E and on the 
parallel momentum k\\\ 

A„(E,k||) =-^ImTr f dV„ G+(r„, r„, E, ky), (1) 

On 

where G'’''(r„, r„, i?, ky) is the retarded Greens function at position r„ in a selected 
atomic sphere of layer n, and the trace is taken over the quantum numbers of the 
total angular momentum. Thus, the BSF is ideal to study the surface states, whereas the 
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three-dimensional bulk band structure is projected onto the two-dimensional Brillouin 
zone (BZ) corresponding to the crystal surface [see Fig{T][b)]. 

3. Band bending and shifting of the Dirac cone by charge dopants 

We now turn our attention to the SSs formed on the surface of Bi 2 Se 3 crystal. In 
this section we discuss our results on the BB effect induced by the electric charge 
accumulation and/or inhomogeneous charge distribution close to the surface. We also 
show that the Dirac cone can be shifted in energy due to the effective electric field 
generated by the deposited surface dopants. Our findings are in good agreement with 
the ARPES experiments [6]. 

Making use of the atomic potentials determined by the self-consistent calculations 
described in Sec. [2] we calculated the layer-resolved BSF in the QLs beneath the 
surface. (The band structure calculated for the bulk crystal is presented in appendix 
Appendix A[) Fig. [2] shows the calculated BSF summed over the layers of the individual 



Figure 2. The BSF given in Eq. o and summed over the layers of the (a) IQL, (b) 
2QL, (c) 3QL and (d) 4QL in the interface region plotted along the KTM cross section 
of the two-dimensional BZ. The projection of the bulk bands onto the two-dimensional 
BZ is shown by a colored areas, while the narrow lines correspond to the bands with 
low dispersion in the z direction, including the SSs. The intensity of the SSs rapidly 
decrease in deeper QLs. The results of the SKKR model are compared to the energy 
bands of the effective model ([3]) using parameter set (a) vq sa 3.55 eVA, A w 128 eVA^, 
1/(2 to) = 0, a = 0 and (b) « 1.63 eVA, A* Ri 108 eVAh 1/(2 to*) « 20 eVAh 

a* Ri 78 eVA^ (see the text for details). The BSF was calculated at complex energies 
with small imaginary part of 0.7 meV. 


QLs. The narrow lines correspond to the bands with low dispersion in the direction, 
including the SSs. Within the bulk band gap the dispersion of the protected SSs form a 
Dirac cone, which is anisotropic in the k\\ plane. The Dirac point is located ~ 250 meV 
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below the bulk conduction minimum, in good agreement with the ARPES measurements 
[0]. The SSs penetrates below the surface up to the third QL, where the intensity of 
the SSs eventually vanishes. The signatures of the bulk band structure, however, can 
be observed in all the QLs up to the topmost one. Thus, the SSs spatially overlap with 
the bulk states, hence they can hybridize. Consequently, one can observe an increased 
width of the SS bands in the vicinity of the bulk bands. 

The SSs that are not hybridized with the bulk states, on the other hand, can be 
described by a 2 x 2 effective Hamiltonian proposed by Fu et ah m up to the third 
order of the momentum Kk = {hkx, hky)\ 

k‘^ X 

H(k) = ^ + Vkik^ay - kyax) + - {k% + ki)a^ . (2) 

Here Vk = vq + ak"^, k± = kx i: iky and crx,y,z stand for the Pauli matrices acting in the 
spin space. Parameters m, Uq, a and A are to be determined either from Erst principle 
calculations or from hts to the experimental data. The energy eigenvalues of the electron 
states are given by the expression 

B±(k) = T ± foik-‘+- 3kir , (3) 

where ± labels the bands above/below the Dirac point. The effective mass m introduces 
an asymmetry between the upper and lower side of the Dirac cone that is, indeed, 
signihcant according to the calculated band structure shown in Fig. O In Refernce |1] 
the parameters of Eq. ([3]) were htted to the experimental ARPES data in the upper side 
of the Dirac cone. Focusing on the parameters describing the most pronounced features 
of the dispersion, namely the Fermi velocity Vq and the parameter A responsible for 
the hexagonal warping, the obtained numerical values were Vq ~ 3.55 eVA, A ~ 128 
eVA^, l/(2m) = 0 and a = 0. Fig. [2]^a) compares the band structure of the effective 
model to the SKKR band structure. As one can see, the effective model with the given 
parameters reproduces the calculated band structure of the SSs above the Dirac point 
well. However, below the Dirac point the agreement of the effective and SKKR model 
is weak. To describe the asymmetry between the upper and lower side of the Dirac 
cone we propose another set of parameters, namely Vq ~ 1.63 eVA, A* ~ 108 eVA^, 
l/(2m*) ^ 20 eVA^ and a* ~ 78 eVA^. Using this parameter set we find that both the 
upper and lower side of the Dirac cone obtained by Eq. ([ 3 ]) is close to the dispersion of 
the SSs calculated within the SKKR framework. 

Comparing the band structure shown in Fig. [21(a) to the ones of Figs.[21(b)-(d), one 
can observe the BB in the topmost QL even without any extra charge dopants added to 
the system. Due to the effective electrostatic field induced by the inhomogeneous charge 
distribution close to the surface, the bulk conduction band is repelled upward by about 
100 meV in the first QL compared to the conduction band minimum in the other QLs 
(see Fig. A BB of similar magnitude on the surface of pristine Bi 2 Se 3 crystal was 
also reported in other DFT calculations [33] . We expect that charge dopants deposited 
on the Bi 2 Se 3 surface further modifies the electrostatic configuration of the top layers. 
Indeed, in Hsieh et. al. [6] the surface of the Bi 2 Se 3 sample was dosed by NO 2 molecules 
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and accomplished to shift the Dirac point towards the condnction band. The presence 
of charged imparities in the sample, like electron donor Se vacancies generated dnring 
the sample preparation process or reactive chemical doping [12], also plays a crncial role 
in the electrostatic properties of the snrface. The slow migration of the Se vacancies 
PEI [13] results in an increased concentration of positively charged impurities close 
to the surface. Thus, a time dependent BB was observed in ARPES experiments 
[21131 Elini El in] where the bands at the surface were gradually bent downward. 
In addition, it has been showed both experimentally la and theoretically [23l [23] that 
a potassium (K) layer deposited on the surface of Bi 2 Se 3 crystal triggers similar BB 
effects on the band structure as the Se vacancies. Thus, the electron donor K adatoms 
can be also described by positively charged impurities close to the surface. 

3.1. Charge dopants on the surface of the Bi^Se^ crystal 

In our calculations we simulated the effect of the surface dopants by a planar capacitor 
situated on the surface of the Bi 2 Se 3 crystal, and charged by 6qs per two-dimensional 
unit cell. For example, the presence of the electron donor Se vacancies or K adatoms 
can be described by a positively charged capacitor. On the other hand, the presence 
of the electron acceptor NO 2 molecules used in the ARPES experiments of Ref. [6] 
can be modeled by a negatively charged capacitor. This is a good assumption because 
the characteristic size of the probed samples are typically much larger than the lattice 
constant and the surface physics can be studied in terms of average physical quantities, 
such as the average concentration of the charged surface dopants. In our model the 
charge of the capacitor corresponds to the average charge transfer between the dopants 
and the surface, while we neglect the inhomogeneities on the atomic length scale. 
Additionally, since we are not conhned by the supercell model of other Erst principle 
calculations [23l [23], we can now study the shift of the Dirac cone as a function of the 
doping in the entire concentration range. 

The electric held of the charged capacitor is accounted for by means of a term in 
the Poisson equation, that is solved self-consistently within the SKKR code. In the 
vacuum the electric held of the capacitor is canceled due to the attracted (or repelled) 
electrons from (to) the bulk reservoir. Fig. [3] shows the induced excess charge (AQ) 
on the individual layers compared to the undoped system. The results were obtained 
for a capacitor located close to the surface, and being charged to Sq = O.Olge, where 
qe ~ 1.6 X 10“^® is the elementary charge unit. In particular. Fig. [3] shows the excess 
charge obtained for two diherent positions of the capacitor. As expected, the charge 
of the capacitor is compensated by the accumulated electrons within the hrst two QLs, 
since the induced excess charge can be hosted only by the unsaturated SSs. It can also 
be noticed that the induced ehective electric held might generate further charge transfer 
between the layers in the low-lying QLs as well. For example, in Fig.[3](a) one can clearly 
observe oscillations of AQ up to the 5th QLs. Since the oscillations are centered around 
AQ = 0, the net charge of the corresponding QLs remains zero within the numerical 
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Figure 3. Charge excess on the individual layers induced by a capacitor charged with 
5q — O.Olge per two-dimensional unit cell. The interface region consists of six QLs, 
surrounded by the vacuum (V) from the left and by the bulk (B) from the right. The 
solid vertical lines indicate the position of the capacitor in each figure. 


precision of our calculations. Our conclusions are consistent with previous theoretical 
works pSl 123] predicting long range oscillations in the charge transfer. However, our 
results indicate that these oscillations becomes suppressed for charge dopants located 
beneath the surface of the Bi 2 Se 3 crystal. In the following we present our result obtained 
for the planar capacitor located at position corresponding to Fig. |3]^a). 

Figure 01 compares the band structures of the doped and clean systems. For 



(b) 0.6 
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Figure 4. Top panels: The BSF at positive surface doping compared to the clean 
system. The BSF was calculated using Eq. © and summed over the atomic layers in 
the (a) IQL and (b) SQL. Bottom panels: Similar to the top panels but for negative 
surface doping. The projection of the bulk bands onto the two-dimensional BZ is shown 
by colored areas, while the narrow lines correspond to the bands with low dispersion in 
the z direction, including the SSs. The BSF was calculated at complex energies with 
small imaginary part of ^ 0.7 meV. 


positively charged capacitor the accumulation of the electrons close to the surface is 
energetically favorable. The states for these accumulated electrons are provided by the 
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downward shift of the unsaturated SS bands [see Figs. 111(a) and (b)]. Moreover, the BB 
of the bulk bands can also be tuned by a surface doping. For a given concentration of 
the surface dopants, for example, the position of the conduction band at the surface can 
be shifted back to that in the bulk. For a negatively charged capacitor one can follow 
analogous reasoning. In this case the bands are shifted upward as shown in Fig. 111(c) 
and (d). 

In previous numerical studies the effective electric held responsible for the BB was 
controlled by a charge transfer originating from a sheet of adatoms located above the 
Bi 2 Se 3 surface [23l Ell ESI [32] . In these calculations additional SSs were identihed in the 
band structure being attributed to the presence of the given adatoms. In the case of 
potassium for example, even quantum well states have been found within the bulk band 
gap where the Dirac cone is situated [23]. Moreover, due to a conventional insulator 
interface on the top of Bi 2 Se 3 , the protected SSs can even be spatially shifted towards 
the bulk of the Bi 2 Se 3 crystal [32]. Since we employ a different approach to model the 
BB effect, these artifacts are entirely absent from our results. As a result, our method, 
enables us to study the evolution of the Dirac cone as a function of the surface doping in 
the entire concentration range, especially in the low concentration limit which is hardly 
accessible for approaches based on a supercell method due to numerical limitations. 

The energy shift of the Dirac cone was successfully demonstrated by ARPES 
measurements as well PET]. As it is shown in Fig. |H the surface bands are shifted 
upward or downward equally in each QL. This behavior is in line with the localized 
nature of the surface states showing no momentum dispersion along the direction 
perpendicular to the surface. On the other hand, the bend of the bulk bands varies 
continuously with the distance measured from the surface. Indeed, in Fig. El(a)-(d) one 
can observe a gradual decrease of the conduction band minimum in the successive QLs. 
The layer resolved BSF (not shown in the manuscript) also conhrms that the energy 
shift of the conduction band minimum varies smoothly from layer-to-layer. 

We now examine the relation between the energy shift of the Dirac cone and the 
concentration of the charge dopants. In case the Fermi energy is located inside the 
bulk band gap the excess charge induced by the surface doping can be hosted only by 
the unsaturated SSs, forming a Dirac cone at the center of the two-dimensional BZ. 
One can then expect, that the energy shift of the Dirac cone is in correspondence with 
the increment or reduction of the occupied electron states required to host the induced 
excess charge. The total charge per two-dimensional unit cell (D) corresponding to the 
electron states on the Dirac cone between energies Ei and E 2 can be calculated as 



(4) 


where 



( 5 ) 
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is the density of states and hvgfk) = |grad]^ii^±(k)| is the gronp velocity. The integral 
path is taken over the constant energy contonr Tg of the SS spectrnm ([3]) at energy 
E. Figure [5l shows the calculated energy shift of the Dirac cone as a function of the 



Figure 5. Energy shift of the Dirac cone predicted by Eq. o as a function of the 
surface doping per a two-dimensional unit cell (solid line). Red crosses represent the 
results obtained by the SKKR calculations. 

deposited excess charge. Surprisingly, we have found a significant difference between 
the predictions of the effective model and the SKKR results. The effective model highly 
overestimates the energy shift of the Dirac cone compared to the results obtained by the 
SKKR model. One can explain this inconsistency by the following arguments. Within 
the effective model we assumed that the surface bands are shifted without any notable 
changes of the dispersion. However, this approximation is valid only at energies close 
to the Dirac point. Indeed, in Fig. 0] the shape of the Dirac cones are similar in the 
doped and undoped cases. One can, however, observe remarkable changes in the surface 
band structure around the crossing points with the other bands when surface dopants 
are present in the system (see Fig. 0]). Consequently, the charge density distribution 
on these parts of the band structure undergo to a signihcant change as well. Thus, the 
low-energy segments of the unsaturated surface bands also play an important role in 
the screening of the surface dopants. 

Comparing our results to other DFT based calculations [231 EU, generally we 
found a good qualitative agreement. The quantitative deviations can be related to 
the differences in the applied physical models. While most of the DFT calculations use 
a supercell approach, our method relies on using semi-infinite bulk and vacuum regions. 
Secondly, in Ref. [23] the charge transfer between the deposited potassium atomic layer 
(K) and the Bi 2 Se 3 surface was controlled by the K-Bi 2 Se 3 distance instead of the K 
concentration on the surface. In this work we found the position of the Dirac point to 
be more sensitive to the amount of charge transfer, resulting in a Dirac cone shift-charge 
transfer relation that is closer to the prediction of the effective model. Our results, on the 
other hand, indicate that the energy shift of the Dirac cone cannot be estimated within 
the effective model given by Hamiltonian (I2|), since a significant portion of the induced 
excess charge is hosted by the low-energy segments of the surface bands. Thus, the 
position of the Dirac point is very much influenced by the treatment of the electrostatic 
potential used in the specific surface calculations. 
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3.2. Spatially distributed charge dopants 

Besides a planar capacitor we also considered a scenario of spatially distribnted charge 
dopants below the snrface, that might be closer to the realistic case of an exposed snrface. 
One can expect that the concentration of the snrface dopants vanish exponentially with 
the distance measnred from the snrface of the crystal. Thus, in our calculations we 
described the concentration of the charged impurities by an exponential function 

Sq^{z) = ^Exp (-Sj . (6) 

The parameter f describes the penetration depth of the dopants into the crystal. The 
case of the planar capacitor can be recovered by the limit ^ ^ 0. In this section we 
discuss the results of our SKKR calculations with finite f parameters. In particular the 
layer dependent excess charge in Fig. [6] was calculated for charge dopants concentrated 
into the hrst two QLs below the surface. As one can see in Fig. [6], the oscillations in the 



Figure 6. Charge excess on the individual layers induced by dopants of total charge 
6q = —O.Olqe and being spatially distributed in the layers close to the surface. The 
smooth blue line describes the concentration of the charge dopants as a function of the 
distance measured from the surface calculated by Eq. (jb)) using ^ ss 5 A. The interface 
region consists of six QLs, surrounded by the vacuum (V) from the left and by the 
bulk (B) from the right. 

excess charge decay on a length scale of 4 — 6 QLs, which is longer that we have seen 
for a planar capacitor located inside the crystal. Our calculations also indicate that the 
decaying length of the oscillations increases with the parameter f, which is consistent 
with our expectations. 

Figured compares the band structures of the spatially doped and clean systems. 
We found very similar results compared to the case when the charge dopants were 
modeled by a planar capacitor on the surface (see Fig. 01). The downward shift of the 
Dirac cone induced by positive charge dopants [see Fig. [71(a) and (b)] is smaller, but 
very close to that induced by the planar capacitor used in the previous section [see Fig. 
111(a) and (b)]. For negative charge dopants, on the other hand, the upward shift of the 
Dirac cone is about the twice of the shift induced by the planar capacitor. However, this 
energy shift is still much less than the predicted value by the effective continuous model. 
Thus, the qualitative conclusions made in the previous section are also applicable for 
the case of spatially distributed charge dopants. 
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Figure 7. Top panels: The BSF at positive spatial doping compared to the clean 
system. The BSF was calculated using Eq. (HD and summed over the atomic layers of 
(a) the IQL and (b) SQL. Bottom panels: Similar to the top panels but for negative 
spatial doping. The projection of the bulk bands onto the two-dimensional BZ is shown 
by a colored areas, while the narrow lines correspond to the bands with low dispersion 
in the z direction, including the SSs. The BSF was calculated at complex energies 
with small imaginary part of ~ 0.7 meV, and ^ r; 5 A. 


4. Summary 

In summary, we have calculated the band structure of Bi 2 Se 3 topological insulator by 
using the SKKR method. In order to examine the effect of the charged impurities on 
the properties of the Bi 2 Se 3 surface, we also calculated the surface band structure in 
the presence of a charged planar capacitor situated close to the surface and for spatially 
distributed dopants under the surface. We have found, that for a Fermi energy located in 
the bulk band gap the induced excess charge is hosted by the unsaturated SSs. Thus, the 
charge of the surface dopants is also screened within the hrst few QLs below the surface. 
In addition, due to the excess charge and the inhomogeneous charge distribution close 
to the surface, the bulk bands undergo a BB effect even in the pristine Bi 2 Se 3 crystal. 
Consequently, the bulk bands becomes bent in the atomic layers close to the surface, 
but one can already recover the properties of the bulk crystal starting from the third 
QL. Our results also indicate, that the BB prohle can be tuned via the deposited surface 
dopants. 

In contrast to the bulk bands, the Dirac cone (being formed by the SS bands inside 
the bulk band gap) becomes shifted in energy due to the deposited surface dopants. 
The magnitude and the direction of this energy shift depends on the concentration and 
on the sign of the deposited dopants. However, this effect cannot be described within 
the scope of the effective continuous model of the SSs. Our self-consistent numerical 
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results showed that besides the Dirac cone the low-energy segments of the surface bands 
also play an important role in the electrostatic properties of the surface. We also found 
that in agreement with recent theoretical studies [23l |23], the charged impurities on 
the surface induce oscillations in the charge density extending deep into the crystal. 
However, our results indicate that in the case when the dopants are located beneath the 
surface, these oscillations became highly suppressed. 

In order to check experimentally our hndings, one needs to independently measure 
the doping concentration on the surface and the energy shift of the Dirac cone. 
We believe that the combination of the scanning tunneling microscope and ARPES 
techniques can serve this purpose. Moreover, the surface doping of the 2GTIs can be 
used to cancel the BB effect on the surface, which is essential to take an advantage of 
the protected SSs in transport devices. Finally, the possibility of tuning the position of 
the Dirac cone might also be of great importance for future experimental applications 
of these materials. 
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Appendix A. The bnlk band structnre of B^Sea componnd 

In this section we determine the band structure of the bulk Bi 2 Se 3 crystal. To this end 
we calculate the layer resolved BSF given by Eq. ([I]) for the bulk system at complex 
energies with a small imaginary part, and summed the contributions of the individual 
atomic and empty layers within one QL. The result of this procedure can be interpreted 



-1 -0.5 0 0.5 

MA-1] 


Figure Al. Projection of the BSF onto the two-dimensional BZ summed over the 
layers of one QL inside the bulk plotted along the KTM cross section of the two- 
dimensional BZ. The blue and red colors represents areas of low and high density of 
states, respectively. The Fermi level is at the top of the valence band. The BSF was 
calculated at complex energies with small imaginary part of ~ 0.7 meV. 

as the projection of the bulk band structure onto the two-dimensional BZ corresponding 
to the crystal intersection formed perpendicularly to the 2 ; direction. (For details see 
the main text.) The calculated BSF along the KTM cross section of the BZ is shown 
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in Fig JAll The high density areas (forming narrow lines close to the extremal points of 
the band structure) indicate electron states of low dispersion along the ^ direction. The 
Fermi level is at the top of the valence band. 

It should be noted that the Fermi level is sensitive to the presence of impurities, 
even at small concentration. Several ARPES measurements indicated that due to the 
bulk dopants the Fermi energy was tuned into the bulk conduction band.[l] According 
to our calculations, the direct bulk band gap at the F point is ~ 447 meV. Our results 
are in good agreement with other theoretical calculations. mEsiiig 

References 

[1] H. Zhang, C. X. Liu, X. L. Qi, X. Dai, Z. Fang and S. C. Zhang, Nature Phys. 5, 438 (2009). 

[2] H.-J. Noh, J. Koh, S.-J. Oh, J.-H. Park, H.-D. Kim, J.D. Rameau, T. Valla, T.E. Kidd, P.D. 

Johnson, Y. Hu, and Q. Li, EPL 81, 57006 (2008). 

[3] D. Hsieh, Y. Xia, D. Qian, L. Wray, F. Meier, J. H. Dil, J. Osterwalder, L. Patthey, A. V. Fedorov, 

H. Lin, A. Bansil, D. Grauer, Y. S. Hor, R. J. Cava, and M. Z. Hasan, Phys. Rev. Lett. 103, 
146401 (2009). 

[4] K. Kuroda, M. Arita, K. Miyamoto, M. Ye, J. Jiang, A. Kimura, E. E. Krasovskii, E. V. Chulkov, 

H. Iwasawa, T. Okuda, K. Shimada, Y. Ueda, H. Namatame and M. Taniguchi, Phys. Rev. Lett. 

105, 076802 (2010). 

[5] Y. Xia, D. Qian, D. Hsieh, L. Wray, A. Pal, H. Lin, A. Bansil, D. Grauer, Y.S. Hor, R.J. Cava, 

M.Z. Hasan, Nat. Phys. 5, 398 (2009). 

[6] D. Hsieh, Y. Xia, D. Qian, L. Wray, J. H. Dil, E. Meier, J. Osterwalder, L. Patthey, 

J. G. Checkelsky, N. P. Ong, A. V. Fedorov, H. Lin, A. Bansil, D. Grauer, Y. S. Hor, R. J. Cava 
and M. Z. Hasan, Nature 460, 1101 (2009). 

[7] Y.H. Wang, D. Hsieh, D. Pilon, L. Fu, D.R. Gardner, Y. S. Lee and N. Gedik, Phys. Rev. Lett. 

107, 207602 (2011). 

[8] A. Herdt, L. Plucinski, G. Bihlmayer, G. Mussler, S. Doring, J. Krumrain, D. Griitzmacher, 

S. Bliigel, C.M. Schneider, Phys. Rev. B 87, 035127 (2013). 

[9] C. Mann, D. West, 1. Miotkowski, Y.P. Chen, S. Zhang, C.-K. Shih, Nature Communications 4, 

2277 (2013). 

[10] J. Honolka, A. A. Khajetoorians, V. Sessi, T.O. Wehling, S. Stepanow, J.-L. Mi, B.B. Iversen, 

T. Schlenk, J. Wiebe, N.B. Brookes, A.I. Lichtenstein, Ph. Hofmann, K. Kern, and 
R. Wiesendanger, Phys. Rev. Lett. 108, 256811 (2012). 

[11] L. Fu, Phys. Rev. Lett. 103, 266801 (2009). 

[12] H.M. Benia, C. Lin, K. Kern, and C.R. Ast, Phys. Rev. Lett. 107, 177602 (2011). 

[13] S.-X. Wang, P. Zhang, S.-S. Li, AIP Advances 3, 052105 (2013). 

[14] M. Bianchi, R. C Hatch, D. Guan, T. Planke, J. Mi, B.B. Iversen and P. Hofmann, Semicond. Sci. 

Technol. 27, 124001 (2012). 

[15] D. Hsieh, J.W. Mclver, D.H. Torchinsky, D.R. Gardner, Y.S. Lee, and N. Gedik, Phys. Rev. Lett. 

106, 057401 (2011). 

[16] D. Hsieh, F. Mahmood, J.W. Mclver, D.R. Gardner, Y.S. Lee, and N. Gedik, Phys. Rev. Lett. 

107, 077401 (2011). 

[17] Z.-H. Zhu, G. Levy, B. Ludbrook, C. N. Veenstra, J. A. Rosen, R. Comin, D. Wong, P. Dosanjh, 

A. Ubaldini, P. Syers, N. P. Butch, J. Paglione, 1. S. Elfimov, and A. Damascelli, Phys. Rev. 
Lett. 107, 186405 (2011). 

[18] J.G. Analytis, J.-H. Chu, Y. Chen, F. Corredor, R.D. McDonald, Z.X. Shen, and I.R. Fisher, Phys. 

Rev. B 81, 205407 (2010). 



Band bending at the surface of Bi 2 Se^ studied from first principles 


15 


[19] G.S. Jenkins, D.C. Schmadel, A.B. Sushkov, H.D. Drew, M. Bichler, G. Koblmueller, M. Brahlek, 

N. Bansal, and Seongshik Oh, Phys. Rev. B 87, 155126 (2013). 

[20] Oleg V. Yazyev, Emmanouil Kioupakis, Joel E. Moore, and Steven G. Louie, Physical Review B 

85, 161101(R) (2012). 

[21] LA. Nechaev, E.V. Chulkov, Physical Review B 88, 165135 (2013). 

[22] J. Henk, A. Ernst, S.V. Eremeev, E.V. Ghulkov, l.V. Maznichenko, and 1. Mertig, Phys. Rev. Lett. 

108, 206801 (2012). 

[23] K. Park, G.D. Beule and B. Partoens et ah. New J. Phys. 15, 113031 (2013). 

[24] T. Forster, P. Kruger, and M. Rohlfing, Phys. Rev. B 91, 035313 (2015). 

[25] K. Govaerts, K. Park, C. De Beule, B. Partoens, and D. Lamoen, Phys. Rev. B 90, 155124 (2014). 

[26] A. Pertsova, G.M. Ganali, New Journal of Phys. 16, 063022 (2014). 

[27] D. Galanakis, T.D. Stanescu, Phys. Rev. B 86, 195311 (2012). 

[28] L.B. Abdalla, L. Seixas, T.M. Schmidt, R.H. Miwa, A. Fazzio, Phys. Rev. B 88, 045312 (2013). 

[29] M. H. Berntsen, O. Gotberg, B. M. Wojek, and O. Tjernberg, Phys. Rev. B 88, 195132 (2013). 

[30] V. N. Men’shov, V. V. Tugushev, E. V. Ghulkov, Jetp Lett. 98 603 (2014). 

[31] V. N. Men’shov, V. V. Tugushev, S. V. Eremeev, P.M. Echenique, and E.V. Ghulkov, Phys. Rev. 

B 91, 075307 (2015). 

[32] G. Wu, H. Ghen, Y. Sun, X. Li, P. Gui, G. Franchini, J. Wang, X.-Q. Ghen, and Z. Zhang, Scientific 

Reports 3, 1233 (2013). 

[33] J. Zabloudil, R. Hammerling, L. Szunyogh, and P. Weinberger, Electron Scattering in Solid Matter 

(Springer, Berlin, 2005). 

[34] D. M. Ceperley and B. J. Alder, Phys. Rev. Lett. 45, 566 (1980). 

[35] S.K. Mishra, S. Satpathy, and O. Jepsen, J. Phys.: Condens. Matter 9, 461 (1997). 

[36] M.P. Lopez Sancho, J.M. Lopez Sancho, and J. Rubio, J. Phys F 14, 1205 (1984). 

[37] M.P. Lopez Sancho, J.M. Lopez Sancho, and J. Rubio, J. Phys F 15, 851 (1985). 

[38] R. Zeller, J. Phys.: Condens Matter 17, 5367 (2005). 

[39] R. Zeller, J. Phys.: Condens Matter 20, 035220 (2008). 

[40] H. Ebert, D. Kodderitzsch, J. Minar, Rep. Prog. Phys. 74, 096501 (2011). 

[41] P. Lloyd, Proc. Phys. Soc. London 90, 207 (1967); P. Lloyd and P.V. Smith, Adv. Phys. 21, 69 

(1972). 

[42] P. Larson, S.D. Mahanti, and M.G. Kanatzidis, Phys. Rev. B 61, 8162 (2000). 

[43] B.M. Fregoso, S. Goh, J. Phy.s.: Condens. Matter 27, 422001 (2015). 




